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Abstract. The theory of complex networks and of disordered systems is used to study the stability and 
dynamical properties of a simple model of material flow networks defined on random graphs. In particular 
we address instabilities that are characteristic of flow networks in economic, ecological and biological 
systems. Based on results from random matrix theory, we work out the phase diagram of such systems 
defined on extensively connected random graphs, and study in detail how the choice of control policies 
and the network structure affects stability. We also present results for more complex topologies of the 
underlying graph, focussing on finitely connected Erdos-Reyni graphs, Small- World Networks and Barabasi- 
Albert scale-free networks. Results indicate that variability of input-output matrix elements, and random 
structures of the underlying graph tend to make the system less stable, while fast price dynamics or strong 
responsiveness to stock accumulation promote stability. 

PACS. 64.60.aq (Networks), 64.60.De (Statistical mechanics of model systems), 89.65.Gh (Economics; 
econophysics, financial markets, business and management) 
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1 Introduction 

The goals of economic policy-makers include the promo- 
tion of economic growth and minimising the effects of 
down-turns. Consequently, understanding the causes of 
business cycles or fluctuations is vital to their efforts. Re- 
search into these causes has itself followed a cyclic pattern, 
often peaking soon after a major economic downturn pQ. 

A successful modern theory in this endeavour is the 
real business cycle (RBC) theory [2J, according to which 
productivity shocks (i.e., changes in oil prices, technology 
and management strategy) induce fluctuations in capital 
accumulation, consumption and other economic indica- 
tors. Hence, the supply side of the economy is responsible 
for business cycles. Another major school of thought con- 
cerning business cycles is the New-Keynesian (NK) view 
|3], where consumer and investor pessimism are respon- 
sible for business fluctuations. A significant requirement 
on both RBC and NK views of economic cycles is that 
they must be 'consistent with the micro-foundations of the 
macro-economy' 0], i.e., that the study of global quanti- 
ties needs to be linked to the behaviour of the microscopic 
constituents of the economy under consideration. 

In this regard, statistical mechanics is an invaluable 
toolbox, as its approach rests on deriving laws for the 
macroscopic behaviour of many-body systems from their 
microscopic rules of engagement. Research along these 
lines has led to a theory for complex networks [51617] . 
Herein, the attention has - to a certain degree - focused on 



the interplay between the topology of underlying interac- 
tion graphs and the robustness of systems of interacting 
agents. The broad scope of these studies includes ecologi- 
cal networks [S], metabolic networks [9] as well as studies 
of man-made networks such as the internet [10] . 

In [11] Helbing et. al. introduce a non-linear model 
for material flow between sectors of an economy. In their 
model, each material, or good i is characterised by its in- 
ventory level, which, by virtue of a flux balance assump- 
tion, depends on the utilisation of good i by other sectors 
in the economy and on the rate with which this good is 
consumed by external agents. Consumption rates, in turn, 
are influenced by the price of the goods via a non-linear 
demand function. The investigation of this model in 
is based on a linear stability analysis of the stoichiomet- 
ric sectorial utilisation matrix, using empirical data. The 
study concludes that the model does not require exoge- 
nous shocks to explain economic fluctuations. Moreover, 
the arrangement of economic units in a coupled network 
may lead to (undesired) global oscillations unless suit- 
able countermeasures are taken, i.e. unless control poli- 
cies are implemented which avert global instability and 
fluctuations. Other models of production networks include 
|l^ll3ll4ll5ll6TT7Tr5j . 

The aim of this present paper is to broaden the scope 
of the findings in [llj by considering ensembles of ran- 
dom input-output matrices instead of one single sample 
of empirical data. We focus on different ensembles of ran- 
domly assigned stoichiometric input-output matrices, and 
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study the stability of a model whose dynamical properties 
are related to the ensemble of random matricies. As nor- 
mally done, we map stability onto the eigenvalue spectra 
of the stoichiometric matrices. For networks with an ex- 
tensive number of connections per node, the spectra can 
be found analytically using techniques from random ma- 
trix theory [19] . The stability of flow networks defined on 
finitely connected Erdos-Reyni (ER) graphs, small-word 
networks (SWN) and scale-free Barabasi- Albert (BA) net- 
works is addressed from numerical diagonalization of the 
corresponding stoichiometric matrices. 

The remainder of the paper is organised as follows. In 
Sect. [2] we introduce the model. Sect. [3] then outlines the 
steps involved in characterising the system's eigenvalues 
and hence its stability. This is followed by a detailed de- 
scription of our results in Sect. |H Finally, in Sect. [5] we 
provide concluding remarks and describe possible exten- 
sions for further work. Technical details concerning the 
calculation of the density and support of eigenvalues of 
Gaussian random matrices with extensive connectivity are 
presented in appendix lAl for completeness. 



2 Model Definitions 

2.1 Set-up of dynamical control policies 

In this section we begin by describing the statistical model 
of material flow networks, which is based on the model 
introduced in [llj . One considers a system consisting of 
i = 1, . . . , N units of production. For each unit i, we de- 
notes by Qi (t) > its rate of operation, given in units of 
delivery cycles per unit time. We also associate a unique 
commodity, denoted as good i, with each unit. Each unit 
interacts with other units in the network by (i) produc- 
ing/delivering goods and (ii) receiving/consuming goods 
from other units. The net flow, at time t, of good i through 
unit j is —aijQj(t), where a.y S K is the difference be- 
tween the amount of good i received minus the amount of 
good i produced by unit j, per delivery cycle. The sign of 
the dij is here chosen as in [11]. The matrix A is related 
to the aggregate Leontief input-output matrix [20]. The 



quantity — ayQj(t) may be positive (net production) or 
negative (net consumption). We will set an = — 1 in the 
following, reflecting the assumption that unit i produces 
one unit of good i per production cycle. If we assume that 
a fraction, Yi(t) > of good i is consumed by sinks out- 
side the system, then the stock of good i at time t, i.e., 
Si(t) is subject to the conservation law 



dSjjt) 
dt 



= -5>iiQj(*)-*K*)- 



(i) 



We now focus on features that affect the rate of production 
of good i at time t. As per [TT] the factors are two-fold: 
(i) if the current stock S%(t) exceeds an optimal or equi- 
librium level, Sf e K + , the rate of production is reduced 
and vice versa. In operations management literature, such 
a strategy, which is referred to as 'Constant Work-in Pro- 
cess' [H], ensures that each economic sector (or factory) 



maintains and mitigates its inventory and backlog; (ii) if 
the rate of stock accumulation is growing, dSi(t)/dt > 0, 
this is an independent reason to reduce Qi(t), and vice 
versa. 

In addition to the aforementioned factors, each unit i 
can potentially be subject to a control strategy that en- 
sures if the rate of production exceeds an optimal value 
Qi € R + known a- priori, then Qi(t) will decrease, and vice 
versa. This strategy rests on the assumption that each pro- 
duction unit is subject to a budget constraint, limiting the 
range of production rates at which it can operate. Putting 
all these features together, we obtain, similarly to [TT] 



1 dQtjt) 
Qi(t) dt 



S? 



St(t) 
Hi dSj(t) 
Si(t) dt 



1+7, 



Qi(t) 



1 



(2) 



where z^, 7, and /x, are sensitivity parameters. Note that 
dQi(t)/dt is here assumed to be proportional to Qi(t) (i.e. 
relative changes in production rates are considered), en- 
suring that Qi(t) > 0, if one starts with non- negative ini- 
tial conditions. 

For large economies there is an additional equilibrat- 
ing mechanism, relating to the price of good i at time 
t, Pi(t) > 0. As in [TT] the factors that affect the price 
are taken to be identical to those affecting the production 
rates, and consequently we use 



dPi(t) 1 Pi(t)dQi(t) 



dt 



at Qi(t) dt 



(3) 



The pre-factor 1 /ai relates to the sensitivity of price change 
to the factors of influence. Specifically, 1/aj is the price- 
responsiveness, i.e. low values of a, imply that prices of 
commodities relax quickly to their equilibrium prices, while 
large values of a, correspond to slow relaxation. In addi- 
tion, it is also assumed that Pi(t) affects the consumption 
Yi(t) via a demand function, fi(Pi{t)), which is non-linear. 
As per standard practice, fi(Pi) is a monotonic decreasing 
function of Pt. We will write 



Yi(t) = K° + Ut)}fi(W)), 



(4) 



where £i(t) € K arc Gaussian random fluctuations and 
Yj° € IR + is the equilibrium value of external consumption 
of good i. The demand function is modeled as 



fi(Pi(t)) = max (0, di - ekPifjt) 



(5) 



where di and di are non-negative real numbers. The equi- 
librium price is Pf G M + . 

We will now specify choices for the network structure, 
i.e. the stoichiometric coefficients a^-, and discuss their re- 
lation to the equilibrium values {Q®, Sf, P®}- In addition 
to these parameters, the model is defined by the variables 

fii,ai,ji,di,di}. The {di} and {di} determine the re- 
sponse of external consumption to changes of price. The 
remaining variables {vi 1 /ii, a*, 7i} lay out the control poli- 
cies of the production units, and determine their dynam- 
ical adaptive behaviour. They are hence the key control 
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parameters an economic policy-maker would adjust so as 
to maximise stability, and to minimise systemic fragility 
and undesired fluctuating or oscillatory behaviour. 

2.2 Structure of interaction matrices 

The freedom to choose appropriate units for Qt(t) allows 
us to re-scale the <2jj such that, the equilibrium fixed point 
(FP) solution of Eq. ([1]) is given by 

1? =-£««' (6) 

3 

i.e. we scale all a,j such that = 1 for all i. A simi- 
lar approach was taken in [11]. We also assume that at 
equilibrium the external consumption of goods is homo- 
geneous for all goods, i.e., Yf = 1. This simplification 
allows us to express Eq. © as 

3 

These conditions ensure that the overall flux of goods, in- 
cluding a non-negative outflow {5^}, is balanced, i.e. that 
no intrinsic creation of material occurs in the system (im- 
possibility of the Land of Cockaigne [22l23j ). 

We furthermore assume Sf = 1, i.e., the desired equi- 
librium stock level of good i corresponds to the net outflow 
Y® of good i per unit time. Again, similar assumptions 
pertaining to Qf — Y® = Sf have been made in [IT]. 

While in [llj a specific input-output matrix, constructed 
from real-world data, was considered, we here focus on 
a synthetic stochastic setting, in which matrix elements 
dij are chosen to be random variables drawn from an 
ensemble. They are held fixed during the course of the 
temporal evolution of the {Qi(t), Si(t), Pi(t)}. In the lan- 
guage of disordered systems theory the matrix elements 
are 'quenched' [24] variables. As mentioned earlier, et^ rep- 
resents the efficiency with which good i is utilised to pro- 
duce good j. Changing is akin to structural changes 
in the production mechanism; adopting new technology, 
for example. It is reasonable to assume that such changes 
occur on a time-scale slower than that of our dynami- 
cal degrees of freedom, hence justifying our approach to 
regard the interaction matrices as quenched random vari- 
ables. The paradigm of networked systems with randomly 
chosen interaction graphs and coupling constants will be 
discussed further below. 

2.2.1 Structure of matrix elements 

We are interested in the case where each unit interacts 
with a fraction of the other units. This consideration may 
be formalised by decomposing 

0>ij — C-ij Uij , (8) 

where cy 6 {0, 1} are quenched connectivity coefficients, 
determining the adjacency matrix of the flow network, and 
describes the amount of good i utilised by unit j. 



The constraint in Eq. ([7]) is satisfied by constructing 
the linear combination of random variables Jij . 

First, drawing the Jy from some ensemble, and taking 
into account an — —1, we set for i ^= j 

a ij — c ij I Jij ~ TT77 J ' ^ 



u j j 

The set Mi in Eq. ^ denotes the elements on row i such 
that Cij = 1. Analogous approaches have been taken in 
|25l26j . This approach breaks down, however, if the ele- 
ments on each row of matrix J are identical. By provid- 
ing a large enough variance for the distribution of Jij we 
ensure that this case is avoided. Constrained random ma- 
trices have previously also been analysed in the context 
of glassy relaxation [57]. In this case, however, the formu- 
lation of the row constraint induced further correlations 
between off-diagonal and diagonal matrix elements. Our 
implementation, as discussed in Appendix [S] avoids this 
issue, hence simplifying further analysis. 

Below, we investigate the stability properties of the 
model for the cases of (i) dilute, but extensively connected 
and (ii) finitely connected ER random networks, (iii) net- 
works exhibiting the small world property and (iv) scale- 
free networks. 

2.2.2 Gaussian dilute ensemble 

We assume Cy = Cji, i.e., we consider the underlying 
network (as defined by the adjacency matrix) to be undi- 
rected. Directionality in the resulting material flow is mod- 
elled by allowing the utilisation parameters Jij to be asym- 
metric, i.e by allowing for cases in which Jij ^ Jji, as we 
will detail below. The connectivity coefficients are drawn 
from the following distribution: 

PiPa) = (l - |p)<Wo + (10) 

where c € K + is the average connectivity per production 
unit. 

In what follows we consider the limit of so-called 'ex- 
treme dilution' [28], where c — > oo and N — > oo, while 
the ratio c/N tends to zero, i.e. c/N — > 0. This may be 
achieved by allowing, for example c ~ C(log N). This as- 
sumption has important implications for the structure of 
the adjacency graph. Firstly, each node in the graph will 
be connected to a vanishing fraction of the total number of 
nodes. Secondly, the length of a typical loop is O(logiV) 
[2"9"j . Thus, taking TV — > oo, the probability of finding 
loops of finite length tends to 0. The environment about 
each node is thus locally tree-like. 

The utilisation parameters Jy are also taken to be 
quenched. To allow for a well-defined thermodynamic limit, 

— > oo and c — > oo, the mean and variance of Jij need 
to scale suitably with c. Specifically, 

Jij = ~~/= x ij ? (11) 
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leads to a mathematically non-trivial regime, where we 
choose the x^ to be Gaussian random variables of zero 
mean and a variance of order O(N ) [30]. Specifically, the 
are independent in pairs and for i ^ j and k ^ I and 
have the following moments: 

(x^) = 0, (xijXu) = Sikdji + rSuSjk- (12) 

The parameter r e [—1,1] describes the degree of corre- 
lations, with fully symmetric interactions given by r = 1. 
The eigenvalue spectra of such Gaussian random matrices 
can be computed fully analytically, see for example [31j for 
results regarding fully connected Gaussian ensembles. The 
extension to the dilute, but extensively connected case re- 
specting the constraint of Eq. ([7]) is straightforward, we 
report some steps of the corresponding calculation in Ap- 
pendix [X] Our investigation of the stability properties of 
such flow networks can hence be carried out analytically 
to a large degree. 

2.2.3 Gaussian finitely-connected ensemble 

We consider to be distributed according to Eq. (fT0|) . 
and continue to take Cjj = Cjj. However, the average con- 
nectivity, c scales as O(N ), while we still consider the 
thermodynamic limit N — > oo. The statistics of the 
are again those indicated in Eqs. (fTTj) and (|12p . 

This scaling of c and J to be O(N ) complicates the 
analytical characterisation of the statistics of eigenvalues. 
Recent efforts [32|33I34|35I36|37I38|39| have lead to an im- 
plicit characterisation of eigenvalue densities in terms of 
population dynamical equations, often used in the context 
of the cavity method. In our analysis of finitely connected 
cases, we do not resort to such tools, but evaluate the 
corresponding eigenvalue statistics via explicit numerical 
diagonalization. 

2.2.4 Small world graphs 

Another ensemble we consider is that where the cy de- 
fine a Small- World Network (SWN) |40f4T] . Under this 
paradigm, one starts from a network in which units are ar- 
ranged on a one-dimensional lattice with periodic bound- 
ary conditions (i.e. a ring). Each unit is then connected to 
2£ {£ 6 N) of its nearest-neighbours (i.e. £ neighbours to 
the right and £ neighbours to the left of the unit on the 
ring). In the context of an economy, 'near', in a stylistic 
sense, models geographic or economic proximity, e.g. two 
production units within a country. The total number of 
undirected links in the system is N I. 

Based on the algorithm proposed in [JD] , starting with 
the first node and its pre-exiting I nearest-neighbours links 
in a clock-wise direction, we re-wire each link with prob- 
ability k, i.e., the nearest-neighbour link is removed and 
replaced by a link to another randomly selected node. In 
an economic context these re-wired links may, for exam- 
ple, represent economic interactions of a given unit with 
units at long 'distances', e.g. in a different country. This 



procedure is iterated for each node. At the end, the to- 
tal number of links is still the same, while the number of 
re- wired, or long-ranged links is kN £. 

An alternative algorithm is that proposed in [41] , where 
starting with the first node, for each of its pre-existing 2£ 
neighbour-interactions, we add an additional link to an- 
other node with probability k' . One then proceeds with 
the second node and so on. At the end of the procedure, 
the expected coordination number per node is 2^(1 + k'). 

Couplings strengths are given by = Jxij, where 
J e R + and the moments of Xij is given by Eq. (fT2"]) . 

Once again, since £ and J scale as O(N ), an analyti- 
cal characterisation of the statistics of eigenvalues of large 
SWN is difficult. We here limit ourselves to numerical di- 
agonalization when addressing small world networks. We 
also compare results obtained for the two construction al- 
gorithms. 

2.2.5 Scale-free networks 

As a final example we study the model on a scale-free net- 
work. To this end we employ a growth process as proposed 
in [32] and construct the underlying adjacency matrix 
as follows: the seed of the growth process is a network 
composed of two nodes, i — 1,2, with c\2 — C21 = 1. At 
each time-step t = 3, . . . , N one further node is added to 
the network, and connects to one of the already existing 
nodes (i = l,...,i — 1) by preferential attachment, i.e. 
the probability of attaching to node i £ {1, . . . , t — 1} is 
proportional to the degree of node i. As shown in [42] this 
leads to a scale- free degree distribution p(k) ~ k~ 3 asymp- 
totically, i.e. in the limit of infinite network size. In our 
simulations this scaling is reproduced faithfully, yielding 
e.g. exponents of —2.9 at system sizes of N = 1000. In our 
analysis below we will use smaller networks of typically 
A*" = 100 nodes for computational reasons (the stability 
analysis entails diagonalization of matrices of size 3N x 3N 
which can be costly if a large number of samples needs to 
be considered). For such sizes a scale-free degree distri- 
bution with a slightly smaller scaling exponent is found. 
We once more take Jy = Jx^, where J e R + and the 
moments of x^ are given by Eq. (|12p . 

2.3 Paradigm of random network models 

The model as we use it here assumes that the interac- 
tions between units in the system constitute a random 
graph in which the presence or absence and the weight on 
each link existing link are quenched random variables, i.e., 
drawn from some distribution and then kept fixed in time. 
Such random structures can, at best, be seen as a minimal- 
ist approximation to real-world flow networks which are 
generally not random in their structure, and which can 
emerge from a growth or evolutionary process in which 
e.g. certain production units go 'extinct' (bankrupt) and 
where new units join over time. 

Nevertheless, studying quenched random structures al- 
lows for a meaningful abstraction of real-world phenom- 
ena and analytical tractability of the mathematical model. 
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Such approaches have been used in a variety of different 
contexts such as neural networks |43|44j , economic activity 
[25145] and ecology [46147) . amongst others. In ecology in 
particular an ongoing debate on the effects of complexity 
on the stability or otherwise has been sparked by the study 
of random community models, and such model systems are 
under active investigation e.g. in [481149] . The random en- 
sembles of graph structures and distributions of couplings 
we use in our work are characterised by parameters such as 
the mean connectivity or variance of the randomly drawn 
elements in the Leontief matrix. The analysis thus allows 
for a specific characterisation of the effects of such param- 
eters on the stability or otherwise of the system, and on 
its dynamical behaviour. In subsequent work one can then 
build on this approach and add more realism by allowing 
the graph itself evolve in time [50151] . 

A further drawback of the present model is the as- 
sumption that equilibrium values Q®,S®, P® are controlled 
externally (e.g. set to unity), and are not outcomes of the 
dynamics itself. Nevertheless, it was shown in [45] that 
the correlation of fixed-points values of the microscopic 
variables with coupling matrix elements in models with 
random interactions may often be ignored for the consid- 
eration of stability properties of random coupling models. 
Our lines of reasoning follow this approach. 



3 Model Solutions 

Here we investigate the properties and solutions to the 
model presented in Sect. [2l As in TT] we henceforth as- 
sume homogeneous model parameters, z/j = v, fa = ju, 
oti = a and 7i = 7 and use the demand function, 



fi(Pi(t)) = f(Pi(t)) = max 0, d - dP % {t) 



(13) 



At the FP we take P° = 1 and assume that fi{P?) = 1. 



3.1 Linear Stability Analysis 

We now investigate how the stability of the fixed-points 
depends on model parameters. To this end we perform a 
linear stability analysis, i.e. the eigenvalues of the Jaco- 
bian of the systems are used to characterise the dynamical 
behaviour of the system when subjected to external per- 
turbations. 

The linearization of the dynamical system about the 
fixed-points is given by 



dSi N 

= - 2J a i3 qj(t) + dpi(t) . 

i—l 



dt 



= — v s 



i(t) 



, \ dsi 
+ 



-TT = VSi(t) + 7ft (t) + LI — 

at a \ at 



(14) 

(15) 
(16) 



We note that the control policy d has dropped out 
of our equations. While the system has dimension 37V, 



its rank is only 2TV. Consequently, TV of its eigenvalues 
vanish. We denote the remaining 27V eigenvalues by Xi t ± 
with i = 1, . . . , TV. Labelling the eigenvalues of the TV x TV 
stoichiometric matrix A by Ej G C, for i = 1 , . . . , TV, we 
obtain, similar to [11] 



Ai ± J A 2 -45, 



where, 



Bi = v 




(17) 

(18) 
(19) 



Due to the TV eigenvalues at zero, the FPs of our sys- 
tem can either be marginally stable or unstable. Such zero 
modes are not unusual in storage systems of the type we 
are considering here, and reflect the effects of an 'inte- 
grating' behaviour of the buffers, which might cause the 
operating point to drift in time [52] , Following [11] we 
hence characterise the stability or otherwise of the system 
in terms of the remaining 27V eigenvalues. Eqs. (T7|) - ([19]) 
relate the eigenvalues A±j of the full 37V x 37V system to 
those of the TV x TV interaction matrix, A. In what fol- 
lows we show that the stability of the full system depends 
on the statistics of the E* only through the support of its 
spectral density. 



3.2 Density of eigenvalue for the Gaussian dilute 
ensemble 

We here establish the average density of eigenvalues, Ej, 
i = 1, . . . , TV, for an ensemble of TV x TV dilute real Gaus- 
sian random matrices, A, defined via Eq. ([9]), where the 
Cy's and Jy's are drawn according to Eq. (fT0|) and Eqs. 
(fTT|l - p^|) . respectively. The density of eigenvalues Ej is 
given by 



/ 1 N 



(20) 



where the (...) denotes an ensemble average. Following 
the lines of Sommers et al [31] , and referring to the real 
and imaginary parts of E by x and y, respectively (i.e. 
E = x + iy), we obtain 

(tt a by 1 , if ((x + 1) I of + (y/b) 2 < J 2 
/'(E) = { (21) 
0, otherwise , 

in the thermodynamic limit, TV — > 00. We have here writ- 
ten a = 1 + r and b = 1 — T. The eigenvalues are 
hence uniformly distributed in the ellipse with major and 
minor axis given by a and &, respectively. These results 
accurately match numerical results, and while we will not 
enter the details of the derivation of Eq. (21) in the main 
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text, Appendix lAl contains some intermediate steps of the 
computation. In particular, with the scaling of couplings 
as chosen above, the result is independent of the connec- 
tivity parameter c. 



3.3 Mapping the eigenvalue support 

The stability or otherwise, and dynamical behaviour of 
the flow system is characterised by the eigenvalue, Eq. 
(jTTJ) , of the 3N x 3N system with the largest real part. 
We denoted this eigenvalue by A m . Working in the ther- 
modynamic limit and considering the map A i— > X±j, 
defined by Eq. (fT7|) . A m is found to lie on the image of 
the boundary of the ellipse defined by Eq. (|2Tj) . Thus in 
order to determine the long-term dynamical behaviour of 
the system only the image of this boundary needs to be 
considered. 

Fig. [1] verifies the validity of this mapping and com- 
pares the analytically obtained boundary of the spectrum 
of the 3N x 3N system against results from direct nu- 
merical diagonalization. The crosses in the figure are from 
numerical diagonalization of the 3iV x 37V system, while 
the solid line is from mapping of the boundary via 

1 



ReAi,, = - [-ReAi ± y/\ A| cos fa/2)) , (22) 
1 



ImA±, t = - (-ImAi ± VI A| sin (<^/2) J , (23) 

where A = Af ~ ABi and ifii — arctan (Im A/Re A)- 
Fig.[T]demonstrates that the analytical theory captures 
the boundary of the spectrum faithfully, and allows one to 
make statements regarding A m . The identification of this 
eigenvalue may be unique only up to complex-conjugation. 

We note that the density of A±,, within the predicted 
support is not uniform. This is due to a non-trivial Jaco- 
bian of the transformation 



p(ReA±, ImAj 



p(x, y) 



d(x, y) 



(9(ReA±, ImA±) 



(24) 



However, randomly sampling Ei according to Eq. (|2"Tj) and 
applying Eq. (|22l) - (l23|) . yields, in the large N limit, a dense 
scattering of A± within the mapped boundary. 

The sequence of spectra shown in Fig. [T] reveals two 
different transitions of the dynamical behaviour of the 
system as the model parameter J, i.e. the variability of 
elements in the coupling matrix, is increased. At small J 
(see panel (a)) one finds ReA m < and ImA m ^ 0, in- 
dicating damped oscillations. As J is increased (see e.g. 
panel (b)) the real part of X m becomes positive with the 
imaginary part still remaining non-zero. This corresponds 
to growing oscillations. As J is increased further A m con- 
tinues to display a positive real part, but it's imaginary 
part vanishes, i.e., panel (d). Hence the system is in an 
exponentially growing phase, where no oscillations are to 
be expected. In order to verify that these eigenvalue dis- 
tributions capture the stability properties and dynamical 
behaviour of the system correctly, we have integrated the 
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Fig. 1. Eigenvalue scatter plots, wherein the crosses are results 
from numerical diagonalization of Eqs. (|14[) -(116 [ ) and the black 
curve is from mapping eigenvalue support Eq. (|21|l via Eqs. 
(I22[) - (|23[ ). For all plots we choose model parameters r = 0.5, 
a = 1, fl = 0.07, v = 1.0, 7 = 0.0 and d = 10.0. In panel 
(a) J = 2.0, (b) J = 5.0, (c) J = 7.5 and (d) J = 8.0. 
System size is N = 400. 



linearised dynamics, Eqs. (|14l I15|16|) numerically, initial- 
ising the system close to its fixed point. Fig. [5] shows the 
resulting behaviour of g = N^ 1 X)i — ^i) as a func- 
tion of time, and results confirm the transitions predicted 
in by the eigenvalue distributions shown in Fig. [T] 

We discuss the phase behaviour in more detail in the 
next section, and focus on studying how the different model 
parameters affect the stability or otherwise of the model. 



4 Results 

In the space of control policies, i.e., a, u, fi, 7, d, and pa- 
rameters r and J characterising the statistics of the input- 
output matrix as well as in dependence on the structure 
of the underlying network, we ask, what policies promote 
stability? Due to the large number of model parameters, 
our investigations necessarily focus on a few specific cuts 
in parameter space. While this is cannot be an exhaustive 
enumeration of effects of all different model parameters, 
we find that the behaviour exhibited in these phases is rich 
and informative regarding the impact of policy changes. 

For a given set of parameters, the system's FP is meta- 
stable if Re A m < or is otherwise unstable. The trajec- 
tories to the FP are damped if ImA TO = 0. If A m has a 
non-zero imaginary part, then the trajectories are charac- 
terised by oscillations. 

We adopt the following notation to distinguish the dif- 
ferent phases: (i) OD: oscillatory decay (ReA r „ < and 
ImA m 7^ 0), (ii) OG: oscillatory growth (ReA m > and 
ImA m 7^ 0), (iii) ED: exponential decay (ReA m < and 
ImA m = 0) and (iv) EG: exponential growth (Re A m > 
and Im A m = 0). 

Finally, for simplicity, we take 7 = 0. In all tested 
cases for 7 > one finds that ImA m = 0. Consequently, 
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Fig. 2. Deviation g of the average stock accumulated in the 
system from the fixed point value as a function of time. Each 
panel corresponds to a different value of J: (a) J = 1.0, (b) 
J = 5.0, (c) J = 7.5 and (d) J = 9.0. The additional model 
parameters were set as 7 = 0.0, r — 0.5, a = 1.0, /1 = 0.07, 
v = 1.0 and d — 10.0. We calculated g by numerically inte- 
grating Eqs. (|14|l - (|16l ) for 100 units and taking the discretized 
time step A = 0.001. 



in that case, the regions in the phase plane will be either 
ED or EG. 



4.1 Gaussian dilute ensemble 

4.1.1 Preliminary observations: effects of couplings strength 
and symmetry of interactions 

Fig. [3] plots the resulting phase boundaries in the (r,J) 
plane. Crossing the lower curve from below Re X m switches 
from negative (FP is meta-stable) to positive (FP is unsta- 
ble). The upper curve separates regions with zero and non- 
zero Im A m , respectively, with oscillatory behaviour found 
below and damped trajectories above the line. Hence, the 
phase space is divided into three regions, OD, OG and 
EG, and the transitions reported in Figs. [T] and [2] corre- 
spond to moving along a vertical line upwards in the phase 
diagram, at fixed r = 0.5. 

For a given degree of symmetry between matrix ele- 
ments, i.e., fixed r, increasing the variability between in- 
teractions, i.e., </, pushes the system from a stable phase 
with damped oscillations to an unstable phase. If the ma- 
trix elements are fully symmetric, r = 1 then the unsta- 
ble phase is always characterised by exponential growth. 
However, for intermediate degrees of symmetry, there ex- 
ists a range of J for which one observes growing oscilla- 
tions. The behaviour depicted extends into the negative r 
region. In particular, as r approaches —1 the upper curve 
becomes increasingly steep and diverges for r = —1, 
where the unstable phase is always characterised by grow- 
ing oscillations. 

Our preliminary observations are: (i) increasing the 
variance of couplings makes the system more unstable and 




Fig. 3. Gaussian dilute model: Phase boundaries in the (_T, J) 
plane. The model parameters are: 7 — 0.0, a = 1.0, fi = 0.07, 
v = 1.0 and d = 10.0. 
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Fig. 4. Gaussian dilute model: Phase boundaries in the {T, J) 
plane. To produce all panels, we chose fixed model parameters 
7 = 0.0, n = 0.07, v = 1.0 and d = 10.0. Panels differs in 
values of a: (a) a — 1.0, (b) a — 0.1, (c) a = 0.05 and (d) 
a = 0.03. 

(ii) as couplings becomes more symmetric the intermedi- 
ate OG phase diminishes and is absent in the fully sym- 
metric case. In what follows we show that these findings 
are robust under the variation of model parameters. 



4.1.2 Effects of price relaxation rate 

Fig. 2] plots phase boundaries in the (-T, J) plane, for dif- 
ferent values of the price-responsiveness policy \ja. We 
find that the qualitative structure of the phase diagram in 
Fig. [3] remains unchanged as a is varied, i.e., for different 
price relaxation regimes. 

For a given r, the critical value J c , separating the 
OD and OG phases, decreases as 1/a is lowered. More- 
over, the J curve, as a function of r, switches from a 
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monotonia increasing to a decreasing function. We con- 
clude that quickly evolving prices, i.e., large l/a promotes 
stability. 

Fig. [5] illustrates this further by plotting phase bound- 
aries in the (l/a, J) plane for different values of r. The 
phase space is once again divided into the three distinct 
regions. The J c curve is an increasing function of l/a. As 
we increase r the OG regions diminish, until for r = 1 
the two phase boundaries coincide exactly with each other, 
excluding the OG region from the parameter space. We 
observe a direct transition from OD to EG. As previously 
mentioned, increasing J moves the system from a stable 
to an unstable phase. 

The results discussed so far are derived from the ana- 
lytically known spectra of Gaussian random matrices and 
large connectivity in the limit of infinite system size. Only 
in this limit are analytical results available. In order to as- 
sess the behaviour of finite systems we consider the prob- 
ability that a finite system finds itself in either of the four 
phases as a function of the model parameters. These prob- 
abilities were computed from a numerical diagonalization 
and a subsequent identification of the largest eigenvalue. 

Fig. [5] reports the relative frequency with which each 
phases occurs (phases not reported in the figure are not 
observed in the simulation). The behaviour of the finite 
system follows that predicted by the theory to a good ac- 
curacy. Discrepancies occur in a parameter regime where 
the theory predicts exponential growth. Here, the analyt- 
ically predicted largest eigenvalue is real and positive. In 
finite systems however, largest eigenvalues with a non-zero 
imaginary part may be found due to finite-size fluctua- 
tions. We observe a similar discrepancy in the lower-right 
panel of Fig. [T] where the theory predicts a real-valued 
largest eigenvalue, but explicit diagonalization at finite 
sizes delivers eigenvalues A m with non-vanishing imagi- 
nary values. We attribute the smearing-out of the OG to 




Fig. 6. Gaussian dilute model: probability of finding the sys- 
tem in a given phase versus variability J of interaction coef- 
ficients. The different symbols give the probabilities that the 
system finds itself in a given phase: OD (circles), OG (dia- 
monds), EG (triangles). The remaining model parameters were 
7 = 0.0, a = 0.2, n = 0.07, v = 1.0 and d = 10.0. Panel 
(a) is for _T = 0.0, while panel (b) shows the case F = 1.0. 
Vertical lines show phase boundaries obtained from the theory 
in the limit of infinite system size. The probabilities were com- 
puted from direct diagonalization of 100 interaction matrices 
at N = 100 and mapping eigenvalues via Eqs. (I22[) -(I23 | I. 



EG transition and co-existence of both phases at large J 
in the upper panel of Fig. [6] to this finite-size effect. 



4.1.3 Sensitivity of production rate to stock accumulation 

The parameter \i is a measure for the sensitivity of a unit's 
rate of production Qi to its rate of stock accumulation 
dSi/dt. In order to characterise the effects of this control 
policy, we compute in Fig. [7] the phase diagram in the 
(r, J) plane for different values of /i. For any degree of 
asymmetry in the couplings, we find that increasing the 
sensitivity of the production rate to stock accumulation 
enhances stability - in panels (a)-(c) the area covered by 
the OD phase increases as we increase /i. For /i > 0.7 we 
begin to observe the ED phase as well. As the sensitivity is 
increased further the lines separating ED-OD and OG-EG 
transitions converge until the system experiences a direct 
transition from the ED to EG phase. 

Finally, we further support our findings on the roles 
played by a and /x in promoting stability, by investigat- 
ing the special case of fully symmetric matrix elements, 
r = 1. Here all Ej are real and in this case one can show 
that A m is real and negative whenever the following two 
conditions are satisfied simultaneously: 



a < 



2 J — 1 ' 



and 



v_ (d/a - x m ) 
< 4 ' 



(25) 
(26) 
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Fig. 7. Gaussian dilute model: effects of control policy /i: 
Phase boundaries in the (_T, J) plane. To produce all panels, we 
chose fixed model parameters 7 = 0.0, a = f.0, v = f.O and 
d = 10.0. Panels differ in the value chosen for fi: (a) p = 0.01, 
(b) (j, = 0.07, (c) /i = 0.5 and (d) fj, = 0.7. 



Here x m is the real part of the eigenvalue E of the Leontief 
matrix, giving rise to the relevant eigenvalue A m determin- 
ing the stability of the system. A condition similar to the 
latter is also found in [11]. For further illustration, Fig. 
[8] shows the resulting phase diagram in the (1/a, v/fx 2 ) 
plane. The vertical line is where a c is equal to the right 
hand side of Eq. ([25)) , separating the unstable EG phase 
at large a from the stable ones at smaller values of a. 
If prices evolve sufficiently fast (1/a > l/ot c ), increasing 
H stabilises the system further (at fixed v) by supressing 
oscillations, as seen in Fig. [8j In ,TT| a related phase di- 
agram is shown, based on stochastic modifications of an 
input-output matrix drawn from real world data. 

Evaluating Eq. (|26|) . which defines the phase lines sep- 
arating the OD from the ED phase, requires the knowl- 
edge of function of the model parameters. This 
relation may in general be somewhat intricate, at r = 1 
however one has —2 J — 1 < x m < 2 J — 1. The two dashed 
lines in the figure correspond to Eq. ([26)) evaluated for val- 
ues of x m at the limits of this interval, hence limiting the 
location of the phase transition separating the ED from 
the OD phase. The actual transition point is found to lie 
between these boundaries (solid line). 



4.2 Gaussian finitely-connected ensemble 

While the results presented for the Gaussian ensemble 
with a large connectivity are rich, the assumption that 
production networks display such topology is admittedly 
a very stylized one, chosen because of the available analyt- 
ical results of the spectra of the corresponding couplings 
matrices. More realistic choices might correspond to net- 
works in which each node is connected, on average, to a 
finite number of other nodes, c. In order to study such 
models we next switch our attention to finitely connected 
ER graphs as introduced in Sec 12.2.31 



Fig. 8. Dilute, but extensively connected random graph: Phase 
boundaries in the (1/a, v/fJ- 2 ) plane, for the case F = 1.0. 
Additional model parameters were set as J — 1.0 and d 
1(1.0. 



Results for the density of eigenvalues for such cou- 
plings matrices are so far only available for fully sym- 
metric couplings and have been expressed using approxi- 
mation schemes, such as the single defect approximation 
(SDA) [53] and effective medium approximation (EMA) 
[55] , Other approaches, motivated by the statistical me- 
chanics analysis of spin-glass type systems express the 
density via cavity equations [37] and as the solution to 
population-dynamics equations [36138] . 

For asymmetric random matrices, however, similar re- 
sults are not yet available. It has nevertheless been argued 
54J that for any finitely connected ER random graph the 
tails of the spectra is characterised by Lifshitz tails; hence 
the entire complex plane serves as the support. While most 
eigenvalues are found concentrated around an origin, there 
are outlier eigenvalues. For finite sized systems these out- 
liers may have a significant impact on the dynamical be- 
haviour of the system 

Fig. E] plots for finite systems at an average connectiv- 
ity of c = 1.0 the probability that the system finds itself 
in a particular configuration, as a function of 1/a. For 
small 1/a, the system is unstable and in the EG phase. 
However, as we increase 1/a and allow for faster price 
relaxation dynamics, the system first enters the unstable 
OG phase and later makes a transition to the stable OD 
phase. When tested against c = 4 and c = 8 we found 
that the location of the transition points vary by only 
1/a = ±0.3. In particular, this point is in good agree- 
ment with that for the fully connected system as given in 
Fig. [5] and depicted in Fig. [H] by the vertical dashed lines. 
This is a result of our scaling of the matrix interaction 
term ay as 1 / y/c. 

If we remove this particular scaling of the couplings 
the behaviour of the system becomes dependent on the 
mean connectivity of the underlying network, see Fig. [TO] 
As c increases, the system tends to become more unsta- 
ble. In the case of full asymmetry a transition between 
an oscillatory decaying phase and an oscillatory unstable 
phase is found. If interactions are fully symmetric, then a 
we increase c we system switches from the OD phase to a 
state dominated by the EG phase with high probability. 
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Fig. 9. Gaussian finitely connected model: as a function of 
1/q we plot the probability that a particular configuration is 
achieved. The different symbols denote the different phases: 
OD (circles), OG (diamonds) and EG (triangles). For inter- 
mediate 1/q the system is in the OG phase, but transits into 
the OD phase on increasing 1/q further. The average connec- 
tivity was c = 1.0. Additional model parameters taken were: 
r = 0.5, J = 20.0, (J, = 0.07, v = 1.0 and d = 10.0. For 
each a, the probabilities were evaluated from numerical diag- 
onalization of 50 matrices with N = 100. 




Fig. 11. Small World Network: as a function of the re- wiring 
rate k we plot the probability that the system finds itself 
in a particular phase. The different symbols give the differ- 
ent phases: OD (circles) and OG (diamonds). The nearest- 
neighbour links I — 3. Additional model parameters taken 
were: J = 1.0, fi = 0.07, v = 1.0, T = 0.0, a = 1.0 and 
d = 10.0. For each k, the probabilities were evaluated from 
numerical diagonalization of 50 matrices with TV = 100. 




Fig. 10. Gaussian model on random graphs: probability the 
system finds itself in a particular phase as a function of the 
average connectivity c. The different symbols give the different 
phases: OD (circles), OG (diamonds) and EG (triangles). Panel 
(a) is for r = 0.0 while panel (b) show the case F — 1.0. 
Additional model parameters taken were: J = 5.0, a = 0.2, 
H — 0.07, v = 1.0 and d = 10.0. The probabilities were 
computed from direct diagonalization of 50 matrices with N = 
100. 




Fig. 12. Small World Network: as a function of the rate k 
with which links are added to an initially local network we 
plot the probability that a stable configuration is achieved. 
The stable state is characterised solely by the OD phase and 
the system makes a transition into the unstable OG phase. 
Additional model parameters taken were: J = 1.0, n = 0.07, 
v = 1.0, r = 0.0, a = 1.0 and d = 10.0. For each k', the 
probabilities were evaluated from numerical diagonalization of 
50 matrices with N = 100. 



4.3 Small world networks 

In considering the SWNs ensemble, one introduces two 
new parameters to the system: (i) number of nearest- 
neighbours links, £ in the initial network, and (ii) the 
amount of long-range links, parametrised by model pa- 
rameters k and explained in Sec. 12.2.41 

Under the first variant for the algorithm, (3D]> each 
nearest-neighbour link is re-wired with probability k. In 
Fig. QT] we show that as k tends to 1, i.e., the network 
moves from a structured state to a random one, the system 
becomes more unstable. Initially the system finds itself in 



the OD phase, but as k increases we observe co-existence 
with the OG phase. 

For the second variant to construct the network, |41j . 
Fig. HH plots the probability a stable state, i.e., P s taUe = 
Pod + Ped is achieved as a function of the rate k! with 
which additional links are introduced into the nearest- 
neighbour network. As k' increases, the system becomes 
more unstable. Similarly, stability is reduced as the num- 
ber of nearest-neighbours, £, is increased. Thus, in accor- 
dance with results in Sec l4.2[ it is mostly the total number 
of links in the system that (all other parameters remaining 
the same) controls the stability of the network. 



Kartik Anand, Tobias Galla: Stability and dynamical properties 



4.4 Scale-free networks 

Results for scale-free networks are shown in Fig. [13] for 
a specific choice of parameters. We show the probability 
Pstabie as a function of the coupling strengths J for both 
BA like networks and ER graphs. In both cases the model 
is set up such that the average degree of nodes is < k 
2. As in previous cases, the system size is N — 100, and 
results in a scaling exponent of approximately 2.4 in the 
BA-case. The remaining parameters are as indicated in 
the figure caption. 

The upper panel reveals that for both types of the 
underlying adjacency matrix, the system exhibits a sta- 
ble phase at low variability of elements of the Leontief- 
matrix, i.e. at low values of J virtually all randomly drawn 
instances are found to be stable. Note that randomness 
here refers to both, the underlying graph and the coupling 
strengths along the links. Both systems exhibit a crossover 
to an unstable phase as matrix elements become more di- 
verse. Crucially however the probability of finding a stable 
instance of the flow network is consistently higher in the 
case of BA-like graphs as compared to ER-networks. Near 
approximately J ss 2 this effect can be significant, rais- 
ing the probability of being stable from about 30 per cent 
(ER) to 60 per cent (BA). This observation is further il- 
lustrated in the lower panel of Fig. [131 where we depict 
the real part of the relevant eigenvalue A m . In the stable 
regime the relaxation time, or resilience of the network 
against perturbations is given by r = l/\ReX m \, and since 
|i?eA m | is consistently higher (in the stable phase) in the 
BA case as compared to the ER case we conclude that the 
BA network appears to be more resilient against external 
fluctuations than the ER graph. Similar statements can 
be made in the unstable phase, where ReX m > for both 
types of networks, but where this real-part is consistently 
larger in the ER case compared to the BA network. Hence 
perturbations to flow-networks defined on ER graphs show 
a much larger growth rate as compared to BA graphs. We 
would here like to stress that the observations presented 
in Fig. [12] are only for one specific combination of model 
parameters, and that the above statements are hence valid 
only pending a systematic investigation of other circum- 
stances. Still, the chose example illustrates that the struc- 
ture of the underlying adjacency matrix can be relevant 
and that scale-free degree distributions may potentially 
promote stability in the context of the present model. 



5 Conclusion 

We have studied the stability and dynamical properties of 
a material flow system defined on a variety of random and 
complex network structures. Results from random matrix 
theory have been used to address models in the extremely 
dilute limit with Gaussian couplings, provided coupling 
strengths are scaled appropriately with the mean connec- 
tivity 0- Our theoretical findings are here in very good 



stable 




Fig. 13. Behaviour of the model defined on a BA network 
as compared to ER graphs. Top panel: probability of finding 
a stable instance, i.e., -P g ^ aD i e as a function of variability of 
matrix elements J. The lower panel depicts the expected value 
of the real part of the relevant eigenvalue A m . Inset in the upper 
panel shows the degree distributions of the studied adjacency 
matrices (solid line is p(k) ~ k~ 2A ). Model parameters are 
N = 100, r = 0,a = 0.2, fj, = 0.01, v = l,d = 8. Each data 
point is obtained from sampling 200 realizations of the flow- 
network, and subsequent direct diagonalisation. Mean degree 
is c = (k) — 2 for both considered types of networks. 



agreement with results obtained from direct numerical di- 
agonalization, and indicate complex phase diagrams in 
dependence on the network structure, symmetry or oth- 
erwise of input and output matrices and most crucially 
on parameters characterising the external control policies 
applied to the network. 

The analysis of Gaussian dilute networks reveals the 
following key findings: (i) Increasing the variance between 
matrix couple elements makes the system unstable, (ii) 
quickly evolving prices are conducive to a stable environ- 
ment, (iii) greater sensitivity of the rate of production 
to stock accumulation can also yield greater stability by 
suppressing oscillations. The role of the coupling symme- 
try is more intricate: if prices adjust slowly or the rate 
of production is weakly sensitive to stock accumulation, 
then increasing the symmetry between the underlying net- 
work couplings makes the system more stable. For quickly 
evolving prices or high sensitivity, however, increasing the 
symmetry results in a more unstable system. 

An investigation of finitely connected networks allows 
us to study the question of how the mean degree of nodes 
in the flow network affects its stability. At fixed connec- 
tivity we find that symmetry in interactions and a high 
degree of variability in interaction strengths reduces the 
extent of the stable region in parameter space, as in the 
fully connected model. Increased connectivity adds to the 



1 However, it should be noted that the assumption needed 
to carry through the RMT as described in this paper is c ^> 1 



rather than sparseness. Hence, the analysis applies to networks 
with non-sparse connectivity with c = O(N) as well. 
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variance of coupling strengths, and hence again promotes 
instability. 

Studying flow dynamics on SWNs reveals that an in- 
creased the re-wiring probability makes the system more 
unstable if prices evolve slowly. Thus, networks with a 
regular structure, as opposed to random networks, pro- 
mote stability. Similarly, increasing the number of nearest- 
neighbour links has a marked impact of making the system 
more unstable for slow price dynamics. Finally, we have in- 
teraction graphs of BA type, and our results indicate that 
the underlying scale-free structures may promote stability. 

The route taken in this paper was to study flow sys- 
tems defined on random networks. Randomness here refers 
to both stochasticity in the strengths of interactions (en- 
tries in the Leontief matrix), but also to the presence or 
absence of links between individual units of production. 
The latter type of stochasticity is a common tool in the 
theory of complex networks, while the former has been 
applied in a variety of contexts, e.g. in ecology |46|47j . 
linear economies [55] . evolutionary game theory |56|57j . 
Real- world production networks are of course not random, 
neither in their structure nor in the magnitude of inter- 
unit dependencies. Still, choosing ensembles of random 
networks allows one to unearth general principles that 
determine the stability or otherwise of such models, e.g. 
our study consistently seems to indicate that an increased 
variability of elements in the Leontief input-output ma- 
trix generally induces instability. This leads to the obvi- 
ous task of identifying analogous measures of variability 
in real-world production networks, and to verify whether 
or not such a correlation between complexity and stability 
can be confirmed. 

It is also legitimate to ask what values the various 
model parameters would take in real-world production 
networks? To sensibly answer this question one must anal- 
yse real-world data, which is beyond the scope of the 
present paper. Statistics for J and r can be obtained 
by analysing real- world data such as used in [11]. Posi- 
tive values of r indicate a tendency towards two-cycles 
in the production network (e.g. A uses B and B uses A). 
Such direct cycles are presumably unlikely and we expect 
small values, r w more realistic. Control policies such 
as the price responsiveness are dynamic quantities. This 
makes their direct measurement difficult, as the temporal 
behaviour of real systems would have to be probed. While 
calibrating the model is a necessary future step, the focus 
of the present paper is on the statistical mechanics anal- 
ysis and phase behaviour of the model. The contribution 
of the present paper is a comprehensive analysis, outlin- 
ing the complex interplay between relevant parameters, 
against which real-world scenarios may be placed. 

Further directions of future research and modelling at- 
tempts include the intricate dynamics of evolving produc- 
tion networks, in which the underlying graph is a function 
of time itself, leading to a system in which discrete de- 
grees of freedom (absence or presence of links) interacts 
with continuous ones (e.g. production rates defined on the 
nodes of the network) . Such systems are known as hybrid 
complex systems [58] . Depending on the separation of time 



scales, different types of dynamical behaviour might then 
to be expected, with the freedom of removing or adding 
links potentially helping to stabilise the model. 
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A Random Matrix Theory 

Following lines of reasoning provided in |31j and employing 
an electrostatic analogy, the starting point to evaluate the 
spectra of eigenvalues is the Green's function, 



G(E) 



N ^ 



1 



E,; - E 



(27) 



where E^ denotes the eigenvalues of matrix A. The real 
and imaginary parts of Eq. (|27[) relate to an electric field, 
with charges at points Ej. We can define a potential, 

0(E,E*) = -1/JV log det{(A T - E*) (A - E)} , (28) 

where E* is the complex-conjugate of eigenvalue E. Eq. 
(|2gf satisfies, for E ^ E;, d<f>/dE = G(E). The average 
density of eigenvalues p was shown to be related to <ft av- 
eraged over an ensemble of matrices A, i.e., the disorder 
average, (<fi), via Poisson's equation, 



1 

47T 



V 



(29) 



where, V = 4<9 2 / <9E<9E*. The righthand-side of Eq. ([29]) 
vanishes if G(E) satisfies the Cauchy-Riemann conditions 
and is an analytical function in the complex plane. In 
other words, we can re-interpret p as the measure of non- 
analyticity in G(E). 

As a first step, one must regularise the logarithm in Eq. 
(|29| by introducing a positive infinitesimal e. This ensures 
that the matrix whose determinant we seek is positive 
definite; consequently one may represent the determinant 
[59j as an integral over complex variables, 




exp 



{-e{z\z) 



(30) 



where {z* , z) = ^ z* z { and M = (A T - E*) (A — E). 

The next step is to perform the average ((. . .)) via the 
replica trick [60], In x = lim n ^o {x n ~ 1) / n, the result 
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of which may be solved in the limit N — > oo via the saddle 
point technique. However, as noted in [61] . the integrand 
that one evaluates has its extremum for order-parameters 
that do not depend on the replica index. As a result, ( (f>) 
may be calculated formally by setting n = 1. 

The average is facilitated by linearising terms quadratic 
in J via a complex Hubbard-Stratonovich transformation, 



oN(cf>) 




exp 



{ - z ) 



(31) 



where the disorder terms are restricted in 
D = cxp{i(z*, (A T - E*)y) + i(y*, (A - E) z) } 

and 



J 



1 



k 



•f j i k 



(32) 
(33) 



for off-diagonal elements, while a,, = — 1. The ele- 
ments Uij are a linear combination of Gaussian random 
variables Xij and we verify that up to leading 0(1), 



0. 



(uijUu) = S itk Sj t i + rSijSj^, (34) 



We perform the average of the Cij and Uij in the usual 
manner; for details refer to [30] . 

We note that in Eq. ([31]) if we rotate all Zi and yi, i.e., 
multiply them with A = e l9 that has unit modulus, the 
integral remains unchanged. More precisely, the terms in 
the exponential in Eq. (|31[) are invariant under this trans- 
formation. In our case, we must take terms of the form 

n-^EM) 2 , n-^M) 2 , A^EiO*) 2 . r'E.te) 2 , 

N 1 Y^,i z t V* an d N J2i z i Vi an to be equal to zero. We 
consequently introduce order-parameters 



(35) 
(36) 



Consequently, Eq. ([31]) reduces to 



oN{4>) 



N r 



n 

fe=i 



d 2 z k d 2 y k 



exp 



N 



J 2 uv 



2ir 2vr 
J 2 |( W 2 + K) 2 ) 



e u 



(37) 



i(E* + l)w - i(E + l)w* 



We introduce the order-parameter definitions Eqs. (|35|) 
- (|36|) into Eq. (l37|) via Dirac 6 functions. This allows us 
to perform the integrals over z k and yk using Gaussian 
identities, which yields 



= V{.. .) expjiVl-! + ~ 2 + ~ 3 }} , (38) 



where, 



t2 

J uv 



- i[(E* + l)w + (E 
= — ln[iuiw — iwiw*] 



j 2 ^[ w 2 + K) 2 ] 

IK], 
f- iw* u; , 



(39) 
(40) 
(41) 



and X>(. . .) denotes the integral over variables the order 
parameters and their "hatted" conjugate variables, which 
are a consequence of a Fourier representation of the Dirac 
5 functions. 

The average potential, ( <fi ) is related to the saddle- 
point value of the function & = Si + S 2 + S3. We 
proceed by eliminating the conjugate variable by requir- 
ing stationarity. Next, we introduce r = uv > and 
eliminate u using the stationarity condition dty/du = 0. 



- i[(E* + \)w + (E 
+ 2 + ln[r — w w*] . 



j 2 |[K) 2 + w 2 ; 



l)w* 



(42) 



In the limit e — *• + , we may distinguish between two 
possibilities; ^ may take a maximum at r = or there 
may be an extremum for some r > 0. Considering first 
r = 0, Eq. (|42|) resolves to an analytical function in the 
domain of E; hence, G(E) = and p = 0. 

In the case one obtains an extremum for r > 0, with 
E = x + iy, we obtain saddle point equations 



1 



w w 



w 



w 



J 2 ' 

i2(x - 1) 

j 2 (i + rj ' 
j 2 (i - r) 



(43) 
(44) 
(45) 



Evaluating (<p) at this saddle point and employing Eq. 
(f29| we obtain the average density, p = l/vr(l - T 2 ), 
which is valid in the region r > 0, i.e., 



(x 



b 2 



< j 



(46) 



where a = 1 + T and b = 1 — r. 
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